#!/usr/bin/env python
###############################################################################
# $ID$
#
# 项目:GDAL2Tiles,Google Summer of Code 2007&2008
# 全局地图切片类
# 用途:将栅格转换为TMS分片,创建KML SuperOverlay EPSG:4326,
# 基于Google地图和OpenLayers生成简单的Html查看器
# 作者:Klokan Petr Pridal,Klokan Dot Cz的Klokan
# 网页:http://www.klokan.cz/projects/gdal2tiles/
#
###############################################################################
# 版权所有(C)2008 Klokan Petr Pridal。版权所有。
###############################################################################

"""
Globalmaptiles.py
切片地图服务(TMS)配置文件中定义的全局地图切片
==============================================================
生成Web上使用的全局瓦片所需的函数。
它包含实现坐标转换的类,用于:
    -GlobalMercator(基于epsg:3857)
        范围：[-20037508, -20037508, 20037508, 20037508]
        适用于谷歌地图、雅虎地图、微软地图兼容的瓦片
    -GlobalGeodetic(基于EPSG:4326)
        范围：[-179.999997, -85.051129, 179.999997, 85.051129]
        适用于OpenLayers底图和Google Earth兼容切片
更多信息,请访问:
    Http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification
    Http://wiki.osgeo.org/wiki/WMS_Tiling_Client_Recommendation
    Http://msdn.microsoft.com/en-us/library/bb259689.aspx
    Http://code.google.com/apis/maps/documentation/overlays.html#Google_Maps_Coordinates
由Klokan Petr Pridal于2008-07-03创作。
Google Summer of Code2008,OSGEO GDAL2Tiles项目。
如果您在产品中使用这个类,请将其翻译成另一种语言
或发现它对您的项目有用,请让我知道。
我的电子邮箱:klokan at klokan dot cz。
我想知道它是在哪里使用的。
类是在开源GDAL许可证(www.gdal.org)下提供的。
"""

import math


class GlobalMercator(object):
    """
    3857的TMS

    在球面墨卡托投影中生成瓦片所需的函数,
    EPSG:3857(EPSG:谷歌,谷歌地图全球墨卡托),EPSG:3785,OSGEO:41001。
    这些瓦片与谷歌地图、微软虚拟地球、雅虎地图、
    英国军械测量局OPENspace API,...
    您可以将它们叠加在这些Web地图应用程序的基础地图上。

    像素和瓦片坐标采用TMS表示法(原点[0,0]在左下角)。
    TMS Global墨卡托瓦片需要哪些坐标转换:
         LatLon      <->       Meters      <->     Pixels    <->       Tile     
     WGS84 coordinates   Spherical Mercator  Pixels in pyramid  Tiles in pyramid
         lat/lon            XY in metres     XY pixels Z zoom      XYZ from TMS 
        EPSG:4326           EPSG:3857                      :                   
         .----.              ---------               --                TMS      
        /      \     <->     |       |     <->     /----/    <->      Google    
        \      /             |       |           /--------/          QuadTree   
         -----               ---------         /------------/                   
       KML, public         WebMapService         Web Clients      TileMapService
    在EPSG:3857中地球的坐标范围是多少?
        [-20037508.342789244,-20037508.342789244,20037508.342789244,20037508.342789244]
        常数20037508.342789244来自以米为单位的地球周长,
        4万公里,坐标原点在范围的中间。
        实际上,你可以这样计算这个常量:2*math.pi*6378137/2.0
        $ECHO 180 85|gdalTransform-S_SRS EPSG:4326-t_SRS EPSG:3857
        Abs(纬度)大于85.05112878的极地地区将被剪裁。
    带有EPSG:3857的棱锥体的缩放级别常量(像素/米)是多少?
        整个区域位于由256x256像素瓦片覆盖的金字塔顶部(缩放=0),
        每个较低的变焦级别分辨率总是除以2
        初始分辨率=20037508.342789244*2/256=156543.03392804062
    TMS和Google Maps/QuadTree瓦片命名约定有什么不同?
        瓦片栅格本身是相同的(相等的范围、投影、像素大小),
        相同的栅格瓦片只有不同的标识。
        TMS中的瓦片从左下角的[0,0]开始计数,id为XYZ。
        Google将原点[0,0]放置在左上角,引用为XYZ。
        微软正在通过网站上定义的四叉树名称来引用瓦片:
        Http://msdn2.microsoft.com/en-us/library/bb259689.aspx
    纬度/经度坐标使用的是WGS84基准,对吗?
        是的,我们提到的所有经纬度都应使用WGS84大地基准。
        嗯,像谷歌地图这样的网络客户正在通过球面墨卡托,
        因此实际上球面上的经度/经度坐标被视为他们在WGS84椭球体上。

    来自MSDN文档:
        为了简化计算,我们使用球面投影形式,而不是
        椭球体的形式。由于投影仅用于地图显示,
        而且不是用来显示数字坐标的,我们不需要额外的精度
        椭球体的投影。球面投影使大约
        在Y方向上0.33%的比例扭曲,这在视觉上并不明显。
    如何在EPSG:3857中创建栅格并使用PROJ.4转换坐标?
        您可以使用标准的地理信息系统工具,如gdalwarp、cs2cs或gdalTransform。
        所有工具都支持-t_srs'epsg:3857'。
        对于其他地理信息系统程序,请检查投影的准确定义:
        欲了解更多信息,请访问http://spatialreference.org/ref/user/google-projection/
        相同的投影编码为EPSG:3785。WKT的定义在官方中EPSG数据库。
      Proj4 Text格式:
        +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
        +k=1.0 +units=m +nadgrids=@null +no_defs
      EPGS:3857 人类可读的WKT格式:
         PROJCS["Google Maps Global Mercator",
             GEOGCS["WGS 84",
                 DATUM["WGS_1984",
                     SPHEROID["WGS 84",6378137,298.2572235630016,
                         AUTHORITY["EPSG","7030"]],
                     AUTHORITY["EPSG","6326"]],
                 PRIMEM["Greenwich",0],
                 UNIT["degree",0.0174532925199433],
                 AUTHORITY["EPSG","4326"]],
             PROJECTION["Mercator_1SP"],
             PARAMETER["central_meridian",0],
             PARAMETER["scale_factor",1],
             PARAMETER["false_easting",0],
             PARAMETER["false_northing",0],
             UNIT["metre",1,
                 AUTHORITY["EPSG","9001"]]]
    """

    def __init__(self, tileSize=256):
        "初始化TMS全球墨卡托金字塔"
        # 瓦片大小
        self.tileSize = tileSize
        # 初始分辨率 = 2 * Π * 地球半径 / 瓦片大小
        self.initialResolution = 2 * math.pi * 6378137 / self.tileSize
        # 156543.03392804062表示瓦片大小为256像素
        # 原点移位 = 2 * Π * 地球半径 / 2
        self.originShift = 2 * math.pi * 6378137 / 2.0
        # 20037508.342789244

    def LatLonToMeters(self, lat, lon):
        "4326转3857"
        # x = 经度 * 原点移位 / 180
        mx = lon * self.originShift / 180.0
        # y = 底数为e的log[tan(90+纬度) * Π / 360] / (Π/180)
        my = math.log(math.tan((90 + lat) * math.pi / 360.0)) / \
            (math.pi / 180.0)
        # y = y * 原点移位 / 180
        my = my * self.originShift / 180.0
        return mx, my

    def MetersToLatLon(self, mx, my):
        "3857转4326"
        # 经度 = (x / 原点移位) * 180
        lon = (mx / self.originShift) * 180.0
        # lat = (y / 原点移位) * 180
        lat = (my / self.originShift) * 180.0
        # 纬度 = 180 / Π * (2 * arctan[e^(lat * Π / 180)] - Π /2 )
        lat = 180 / math.pi * \
            (2 * math.atan(math.exp(lat * math.pi / 180.0)) - math.pi / 2.0)
        return lat, lon

    def PixelsToMeters(self, px, py, zoom):
        "像素坐标转3857"

        res = self.Resolution(zoom)
        mx = px * res - self.originShift
        my = py * res - self.originShift
        return mx, my

    def MetersToPixels(self, mx, my, zoom):
        "3857转像素坐标"
        # 指定级别的分辨率
        res = self.Resolution(zoom)
        # 像素x = (x + 原点移位) / 分辨率
        px = (mx + self.originShift) / res
        # 像素y = (y + 原点移位) / 分辨率
        py = (my + self.originShift) / res
        return px, py

    def PixelsToTile(self, px, py):
        "像素坐标转瓦片行列号"
        # 瓦片X = 向上取整[ 像素x / 向下取整(瓦片大小)] -1
        tx = int(math.ceil(px / float(self.tileSize)) - 1)
        # 瓦片y = 向上取整[ 像素y / 向下取整(瓦片大小)] -1
        ty = int(math.ceil(py / float(self.tileSize)) - 1)
        return tx, ty

    def PixelsToRaster(self, px, py, zoom):
        "将像素坐标的原点移动到左上角"

        mapSize = self.tileSize << zoom
        # 像素x = 像素x
        # 像素y = (瓦片大小 * 2^zomm) - py
        return px, mapSize - py

    def MetersToTile(self, mx, my, zoom):
        "3857转瓦片行列号"
        # 3857转像素坐标
        px, py = self.MetersToPixels(mx, my, zoom)
        # 像素坐标转瓦片xy
        return self.PixelsToTile(px, py)

    def TileBounds(self, tx, ty, zoom):
        '''
            瓦片行列号转3857边界
         (minx, miny)   tx       tx+1
                     ty  --------
                        |        |
                        |   256  |
                        |        |
                   ty+1  --------  (maxx, maxy)

        '''
        # 瓦片XY转3857坐标
        minx, miny = self.PixelsToMeters(
            tx*self.tileSize, ty*self.tileSize, zoom)
        maxx, maxy = self.PixelsToMeters(
            (tx+1)*self.tileSize, (ty+1)*self.tileSize, zoom)
        return (minx, miny, maxx, maxy)

    def TileLatLonBounds(self, tx, ty, zoom):
        "瓦片行列号转4326边界"

        bounds = self.TileBounds(tx, ty, zoom)
        # 墨卡托XY转经纬度
        minLat, minLon = self.MetersToLatLon(bounds[0], bounds[1])
        maxLat, maxLon = self.MetersToLatLon(bounds[2], bounds[3])

        return (minLon, minLat, maxLon, maxLat)

    def Resolution(self, zoom):
        "给定变焦级别的分辨率(米/像素)(以赤道测量)"

        # return (2 * math.pi * 6378137) / (self.tileSize * 2**zoom)
        # (2 * Π * 地球半径) / (瓦片大小 * 2^zoom)
        return self.initialResolution / (2**zoom)

    def ZoomForPixelSize(self, pixelSize):
        "像素大小转缩放级别"

        for i in range(30):
            if pixelSize > self.Resolution(i):
                return i-1 if i != 0 else 0  # We don't want to scale up

    def GoogleTile(self, tx, ty, zoom):
        "TMS切片转Google切片"

        # 坐标原点从范围的左下角移动到左上角
        # 瓦片X , (2^zoom -1) - 瓦片y
        return tx, (2**zoom - 1) - ty

    def QuadTree(self, tx, ty, zoom):
        "TMS切片换微软QuadTree切片"

        quadKey = ""
        # 瓦片y = (2^zoom -1) - 瓦片y
        ty = (2**zoom - 1) - ty
        # 从zoom到0 循环递减1 (四叉树)
        for i in range(zoom, 0, -1):
            digit = 0
            mask = 1 << (i-1)
            if (tx & mask) != 0:
                digit += 1
            if (ty & mask) != 0:
                digit += 2
            quadKey += str(digit)

        return quadKey

# ---------------------


class GlobalGeodetic(object):
    """
    4326的TMS切片

        生成板卡投影中的全局块所需的函数, EPSG:4326, "unprojected profile"
        此类切片与Google Earth兼容(与任何其他EPSG:4326栅格一样)
        你可以在OpenLayers底图上叠加这些瓦片。

        像素和瓦片坐标采用TMS表示法(原点[0,0]在左下角)。

    我们需要对TMS全球大地测量切片进行哪些坐标转换?
        全局大地测量切片使用4326大地坐标(纬度、经度)
        直接作为平面坐标XY(也称为未投影或版面Carre)
        我们只需要缩放到像素金字塔和切割到瓦片。
        金字塔顶层有两块瓦片,所以它不是方形的,而是矩形的。
        区域[-180,-90,180,90]缩放到512x256像素。
        TMS的左下角有坐标原点(用于像素和瓦片)。
        栅格采用EPSG:4326格式,因此与Google Earth兼容。
         LatLon      <->      Pixels      <->     Tiles     
     WGS84 coordinates   Pixels in pyramid  Tiles in pyramid
         lat/lon         XY pixels Z zoom      XYZ from TMS 
        EPSG:4326                                           
         .----.                ----                         
        /      \     <->    /--------/    <->      TMS      
        \      /         /--------------/                   
         -----        /--------------------/                
       WMS, KML    Web Clients, Google Earth  TileMapService
    """

    def __init__(self, tileSize=256):
        # 瓦片大小
        self.tileSize = tileSize

    def LatLonToPixels(self, lat, lon, zoom):
        "4326换像素坐标"
        # 分辨率 = 180 / 256 / 2^zoom
        res = 180 / 256.0 / 2**zoom
        # 像素x = (180 + 纬度) / 分辨率
        px = (180 + lat) / res
        # 像素y = (180 + 经度) / 分辨率
        py = (90 + lon) / res
        return px, py

    def PixelsToTile(self, px, py):
        "像素坐标转瓦片行列号"
        # 瓦片x = 向上取整[ 像素x / 向下取整(瓦片大小) ] -1
        tx = int(math.ceil(px / float(self.tileSize)) - 1)
        # 瓦片y = 向上取整[ 像素y / 向下取整(瓦片大小) ] -1
        ty = int(math.ceil(py / float(self.tileSize)) - 1)
        return tx, ty

    def Resolution(self, zoom):
        "给定变焦级别的分辨率(弧度/像素)(以赤道测量)"
        # 分辨率 = 180 / 256 / 2^zoom
        return 180 / 256.0 / 2**zoom
        # return 180 / float( 1 << (8+zoom) )

    def TileBounds(tx, ty, zoom):
        "瓦片行列号转4326边界"
        # 分辨率 = 180 / 256 / 2^zoom
        res = 180 / 256.0 / 2**zoom
        # (
        #  瓦片x * 256 * 分辨率 -180,
        #  瓦片y * 256 * 分辨率 -90,
        #  (瓦片x + 1) * 256 * 分辨率 -180,
        #  (瓦片y + 1) * 256 * 分辨率 -90
        # )
        return (
            tx*256*res - 180,
            ty*256*res - 90,
            (tx+1)*256*res - 180,
            (ty+1)*256*res - 90
        )


if __name__ == "__main__":
    import sys
    import os

    def Usage(s=""):
        print(
            "用法: globalmaptiles.py [-profile 'mercator'|'geodetic'] zoom lon lat [lonmax latmax]")
        print("\n")
        if s:
            print(s)
            print("\n")
        print("此实用程序为给定的WGS84纬度/经度坐标(或边界框)打印切片列表")
        print("覆盖特定区域的。切片在给定的'profile'中(默认是谷歌地图的'mercator')")
        print("而在给定的金字塔中,则是'zoom'")
        print("对于每一块瓦片,都会打印一些信息,包括EPSG:3857和EPSG:4326中的bbox")
        sys.exit(1)

    profile = 'mercator'
    zoomlevel = None
    lat, lon, latmax, lonmax = None, None, None, None
    boundingbox = False

    argv = sys.argv
    # argv = argv + ['-profile', 'geodetic', '1', '40.426', '117.119']
    i = 1
    while i < len(argv):
        arg = argv[i]

        if arg == '-profile':
            i = i + 1
            profile = argv[i]

        elif zoomlevel is None:
            zoomlevel = int(argv[i])
        elif lon is None:
            lon = float(argv[i])
        elif lat is None:
            lat = float(argv[i])
        elif lonmax is None:
            lonmax = float(argv[i])
        elif latmax is None:
            latmax = float(argv[i])
        else:
            Usage("错误:参数太多")
        i = i + 1

    if profile not in ['mercator', 'geodetic']:
        Usage("错误:对不起，给定的 profile 尚未实现")

    if zoomlevel == None or lat == None or lon == None:
        Usage("错误:至少指定'zoom'、'lat'和'lon'")
    if latmax is not None and lonmax is None:
        Usage("")

    if latmax != None and lonmax != None:
        if latmax < lat:
            Usage("错误:必须同时提供'latmax'和'lonmax'")
        if lonmax < lon:
            Usage("错误:'lonmax'必须大于'lon'")
        boundingbox = (lon, lat, lonmax, latmax)
    tz = zoomlevel

    if profile == 'geodetic':
        geodetic = GlobalGeodetic()
        px, py = geodetic.LatLonToPixels(lat, lon, tz)
        tminx, tminy = geodetic.PixelsToTile(px, py)
        if boundingbox:
            mpx, mpy = geodetic.LatLonToPixels(latmax, lonmax, tz)
            tmaxx, tmaxy = geodetic.PixelsToTile(mpx, mpy)
        else:
            tmaxx, tmaxy = tminx, tminy
        if boundingbox:
            print("\tTMS 4326瓦片范围 [z,minX~maxX,minY~maxY]:",
                  "[%s, %s~%s, %s~%s]" % (tz, tminx, tmaxx, tminy, tmaxy))
            wgsboundsMin = geodetic.TileBounds(tminx, tminy, tz)
            wgsboundsMax = geodetic.TileBounds(tmaxx, tmaxy, tz)
            print("\tWGS84瓦片经纬度范围:", wgsboundsMin[:2], wgsboundsMax[-2:])
        else:
            for ty in range(tminy, tmaxy+1):
                for tx in range(tminx, tmaxx+1):
                    tilefilename = "%s/%s/%s" % (tz, tx, ty)
                    print("\tTMS WGS84 (z/x/y):", tilefilename)
    if profile == 'mercator':
        mercator = GlobalMercator()

        mx, my = mercator.LatLonToMeters(lat, lon)
        print("\tESPG:3857坐标lat lon:", mx, my)
        tminx, tminy = mercator.MetersToTile(mx, my, tz)
        if boundingbox:
            mx, my = mercator.LatLonToMeters(latmax, lonmax)
            print("\tESPG:3857的maxlat maxlon: ", mx, my)
            tmaxx, tmaxy = mercator.MetersToTile(mx, my, tz)
        else:
            tmaxx, tmaxy = tminx, tminy
        if boundingbox:
            print("\n\tTMS瓦片范围 [z,minX~maxX,minY~maxY]:",
                  "[%s, %s~%s, %s~%s]" % (tz, tminx, tmaxx, tminy, tmaxy))

            gminx, gmaxy = mercator.GoogleTile(tminx, tminy, tz)
            gmaxx, gminy = mercator.GoogleTile(tmaxx, tmaxy, tz)
            print("\tGoogle瓦片范围 [z,minX~maxX,minY~maxY]:",
                  "[%s, %s~%s, %s~%s]" % (tz, gminx, gmaxx, gminy, gmaxy))
            if tz > 0:
                quadkeyUL = mercator.QuadTree(tminx, tminy, tz)
                print("\t微软 Quadkey:", quadkeyUL, '(', int(quadkeyUL, 4), ')')

            boundsMin = mercator.TileBounds(tminx, tminy, tz)
            boundsMax = mercator.TileBounds(tmaxx, tmaxy, tz)
            print("\n\t3857瓦片经纬度范围: ", boundsMin[:2], boundsMax[-2:])

            wgsboundsMin = mercator.TileLatLonBounds(tminx, tminy, tz)
            wgsboundsMax = mercator.TileLatLonBounds(tmaxx, tmaxy, tz)
            print("\tWGS84瓦片经纬度范围:", wgsboundsMin[:2], wgsboundsMax[-2:])
        else:
            for ty in range(tminy, tmaxy+1):
                for tx in range(tminx, tmaxx+1):
                    print("\tTMS 3857 (z/x/y):", "%s/%s/%s" % (tz, tx, ty))

                    gx, gy = mercator.GoogleTile(tx, ty, tz)
                    print("\tGoogle(Z, X, Y):", "%s/%s/%s" % (tz, gx, gy))

                    if tz > 0:
                        quadkey = mercator.QuadTree(tx, ty, tz)
                        print("\t微软Quadkey:", quadkey)
                    bounds = mercator.TileBounds(tx, ty, tz)
                    print("\n")
                    print("\tEPSG:3857 XY范围: ", bounds)
                    wgsbounds = mercator.TileLatLonBounds(tx, ty, tz)
                    print("\tWGS84经纬度范围:", wgsbounds)
                    print("\tgdalwarp -ts 256 256 -te %s %s %s %s %s %s_%s_%s.tif" % (
                        bounds[0], bounds[1], bounds[2], bounds[3],
                        "<your-raster-file-in-epsg3857.ext>", tz, tx, ty))
                    print
